perm filename COMFM.F4[JC,MUS] blob
sn#075946 filedate 1973-12-05 generic text, type T, neo UTF8
00100 C This is taken from the "Handbook of Mathematical Functions"
00200 C by Abramowitz and Stegun, Dover press S1272, 1965.
00300 C Chapter 9 is on Bessel Functions. The polynomial
00400 C approximations are found in section 9.4, pages 369
00500 C and 370. The recursion formula for generating higher
00600 C orders from J0 and J1 is found at the bottom of
00700 C table 9.1, page 390
00800
00900 SUBROUTINE BESCOEF(MI,X,I)
01000 REAL MI,J0,J1
01100 COMMON FREQ(3,0/200)
01200 IF(I.GT.1)GO TO 30
01300 IF(MI.GE.3.0)GO TO 10
01400 xs=(MI/3)**2
01500 j0=1.0+xs*(-2.2499997+xs*(1.2656208+xs*(-0.3163866+
01600 1xs*(0.0444479+xs*(-0.0039444+xs*0.00021)))))
01700 j1=MI*(.5+xs*(-0.56249985+xs*(0.21093573+
01800 1xs*(-0.03954289+xs*(0.00443319+xs*(-0.00031761+
01900 2xs*0.00001109))))))
02000 GO TO 20
02100 10 xs=3.0/MI
02200 j0=(1.0/sqrt(MI))*(0.79788456+xs*(-0.00000077+
02300 1xs*(-0.0055274+xs*(-0.00009512+xs*(0.00137237+
02400 2xs*(-0.00072805+xs*0.00014476))))))*
02500 3cos(MI-0.78539816+xs*(-0.04166397+
02600 4xs*(-0.00003954+xs*(0.00262573+
02700 5xs*(-0.00054125+xs*(-0.00029333+
02800 6xs*0.00013558))))))
02900 j1=(1.0/sqrt(MI))*(0.79788456+xs*(0.00000156+
03000 1xs*(0.01659667+
03100 1xs*(0.00017105+xs*(-0.00249511+xs*(0.00113653-
03200 2xs*0.00020033))))))*
03300 3cos(MI-2.35619449+xs*(0.12499612+xs*(0.0000565+
03400 4xs*(-0.00637879+xs*(0.00074348+xs*(0.00079824-
03500 5xs*0.00029166))))))
03600 20 FREQ(2,0)=J0
03700 FREQ(2,1)=J1
03800 FREQ(2,2)=J1
03900 RETURN
04000 30 Y=I
04100 IF(MI.LE.0.0000001)GO TO 50
04200 IJ=I*2
04300 X=((2.0*(Y-1.0))/MI)*FREQ(2,IJ-2)-FREQ(2,IJ-4)
04400 RETURN
04500 50 X=0
04600 RETURN
04700 END
00100 COMMON FREQ(3,0/200)
00150 IT=-1.0
00200 ACCEPT 102,CAR,XMOD,CINDEX
00300 MAX=CINDEX+5.0
00400 DO 100 J=0,MAX*2,2
00500 I=J/2
00600 NORDER=I
00700 IF(I.NE.1)CALL BESCOEF(CINDEX,COEF,NORDER)
00800 IF(I.GT.0)GO TO 20
00900 FREQ(1,I)=CAR
01000 FREQ(3,I)=0.0
01100 GO TO 30
01200 20 XI=I
01300 YMOD=XI*XMOD
01400 FREQ(1,J-1)=CAR+YMOD
01500 FREQ(1,J)=CAR-YMOD
01600 IF(I.EQ.1)GO TO 10
01700 FREQ(2,J-1)=COEF
01800 FREQ(2,J)=COEF
01900 10 FREQ(3,J-1)=I
01950 IF(IT)I=-I
02000 FREQ(3,J)=I
02050 IT=-IT
02100 30 CONTINUE
02200 100 CONTINUE
02300 PRINT 101,((FREQ(K,L),K=1,3),L=0,MAX*2)
02400 101 FORMAT(1H 20(3F/))
02500 102 FORMAT(3F)
02600 END